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Abstract. As a model of more general contour integration problems we consider 
the numerical calculation of high-order derivatives of holomorphic functions using 
Cauchy's integral formula. Bornemann (2011) showed that the condition number 
of the Cauchy integral strongly depends on the chosen contour and solved the 
problem of minimizing the condition number for circular contours. In this paper we 
minimize the condition number within the class of rectangular paths of step size h 
using Provan's algorithm for finding a shortest enclosing walk in weighted graphs 
embedded in the plane. Numerical examples show that optimal rectangular paths 
yield small condition numbers even in those cases where circular contours are 
known to be of not much use, such as for functions with branch-cut singularities. 



1. Introduction 

To escape from the ill-conditioning of difference schemes for the numerical 
calculation of high-order derivatives, numerical quadrature applied to Cauchy's 
integral formula has on various occasions been suggested as a remedy (for a survey 
of the literature, see Bornemann 2011). To be specific, we consider a function / 
that is holomorphic on a complex domain D 3 0; Cauchy's formula gives 1 

f{n){0) = lkJ T z ~ n ~ lf{z)dz (1) 

for each cycle T C D that has winding number ind(r;0) = 1. If T is not carefully 
chosen, however, the integrand tend to oscillate at a frequency of order 0(n _1 ) 
with very large amplitude (Bornemann 2011, Fig. 4). Hence, in general, there is 
much cancelation in the evaluation of the integral and ill-conditioning returns 
through the backdoor. The condition number of the integral 2 is (Deuflhard and 
Hohmann 2003, Lemma 9.1) 

and r should be chosen as to make that number as small as possible. Since the 
denominator is, by Cauchy's theorem, independent of T, we have to minimize 



i(r) = jT| z |-»- 1 |/( z )|d| z |. 



(2) 



2010 Mathematics Subject Classification. 65E05, 65D25; 68Rio, 05C38. 

1 Without loss of generality we evaluate derivatives at z = 0. 

2 Given an accurate and stable quadrature method, that condition number actually yields, by 

# loss of significant digits « log 10 ic(T,n), 

an estimate of the error caused by round-off in the last significant digit of the data (i.e., the function /). 
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Figure l. Rectangular path of step size h (taken from Lang 1999, Fig. IV.13). 

Bornemann (2011) considered circular contours of radius r; he found that there is 
a unique r* = r(n) solving the minimization problem and that there are different 
scenarios for the corresponding condition number K*(ft) as n —> 00: 

• k* in) — >■ oo, if f is in the Hardy space H 1 ; 

• limsup n ^ 00 K*(n) ^ M, if / is an entire function of completely regular 
growth which satisfies a non-resonance condition of the zeros and whose 
Phragmen-Lindelof indicator possesses M maxima (a small integer). 

Hence, though those (and similar) results basically solve the problem of choosing 
proper contours for entire functions, much better contours have to be found for the 
class H 1 . Moreover, the restriction to circles lacks an algorithmic flavor that would 
point to more general problems depending on the choice of contours, such as the 
numerical solution of highly-oscillatory Riemann-Hilbert problems (Olver 2011). 
In this paper, we solve the contour optimization problem within the more 
general class of rectangular pathes of step size h (see Fig. 1) as they are known 
from Artin's proof of the general, homological version of Cauchy's integral theorem 
(Lang 1999, IV.3). Such paths are composed from horizontal and vertical edges 
taken from a (bounded) rectangular grid O^ C D of step size h. Now, the weight 
function (2), being additive on the abelian group of path chains, turns the grid O^ 
into an edge-weighted graph such that each optimal rectangular path W* becomes 
a shortest enclosing walk (SEW); "enclosing" because we have to match the winding 
number condition ind(W*;0) = 1. An efficient solution of the SEW problem for 
embedded graphs was found by Provan (1989) and serves as a starting point for 
our work. 

Outline of the Paper. In Section 2 we discuss general embedded graphs in which 
an optimal contour is to be searched for; we discuss the problem of finding a 
shortest enclosing walk and recall Provan's algorithm. In Section 3 we discuss 
some implementation details and optimizations for the problem at hand. Finally, 
in Section 4 we give some numerical examples and compare with the optimal 
circles obtained by Bornemann (2011). 

2. Contour Graphs and Shortest Enclosing Walks 

By generalizing the grid O^, we consider a finite graph G = (V,E) built from 
vertices V C D and edges E that are smooth curves connecting the vertices within 
the domain D. We write uv for the edge connecting the vertices u and v; by (2), its 
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weight is defined as 

d(uv)= / Izl-'-Vtolilzl. (3) 

JUV 

A walk W in the graph G is a closed path built from a sequence of adjacent edges, 
written as (where + denotes joining of paths) 

W = V\V 2 + V 2 V 3 + • • • + VmVi', 

it is called enclosing the obstacle if the winding number is ind(W;0) = 1. The set 
of all possible enclosing walks is denoted by EL As discussed in §1, the condition 
number is optimized by the shortest enclosing walk (not necessarily unique) 

W* = argmind(W) 

wen 

where, with W = V\Vi + v^v^ H Y v m v\ and v m +\ = v\, the total weight is 

m 

d(W)^£d(v j v j+1 ). 

7=1 

The problem of finding such a SEW was solved by Provan (1989): The idea is that 
with V u ,v denoting a shortest path between u and v, any shortest enclosing walk 
W* =w\W2 + W2W3 -f i- w m W\ can be cast in the form (Provan 1989, Thm. 1) 

W* = V Wl/Wj + WjWj +1 + Vw j+1/ w t 

for at least one ;'. Hence, any shortest enclosing walk W* is already specified by 
one of its vertices and one of its edges; therefore 

w* e n = {Vu,v + vw + V w ,u :u e v,vw e e}. 

Provan's algorithm finds W* by, first, building the finite set IT; second, by removing 
all walks from it that do not enclose z = 0; and third, by selecting a walk from 
the remaining candidates that has the lowest total weight. Using Fredman and 
Tarjan's (1987) implementation of Dijkstra's algorithm to compute the shortest 
paths V u ,vf the run time of the algorithm is known to be (Provan 1989, Corollary 2) 

0(|V||E| + |V| 2 log|V|). 

3. Implementation Details 

3.1. Edge Weight Calculation. Using the edge weights d(uv) requires to approxi- 
mate the integral in (3). Since not much accuracy is needed, a simple trapezoidal 
rule with two nodes is generally sufficient: 



d(uv)= j |z|-( n+ %(z)|d|z| 

Juv 

(d{u)+d{v)) =d(uv) 



1 uv 

- v\ 



2 
with the vertex weight 

d( z ) = | z |-(»+i)|/( z )|. (4) 

Although d(uv) will typically have an accuracy of not more than just a few bits, 
we have not encountered a single case, in which a more accurate computation of 
the weights would have resulted in a different SEW W* . 
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3.2. Reducing the size of ft. As described in Section 2, Provan's algorithm starts 
by building a walk for every pair (v, e) £ V x E and then proceeds by selecting the 
best enclosing one. A simple heuristics, which worked well for all our test cases, 
helps to considerably reduce the number of walks to be processed: Let 

v* = argmind(i7) 

vev 

and define W^ as a SEW subject to the constraint 

W^ e fl^ = {V v *,u + uw + V W/V * - uw e E}. 

Obviously W* and W^ do not need to agree in general, as v* does not have to be 
traversed by W* . However, since v* is the vertex with lowest weight, both walks 
differ mainly in a region that has no, or very minor, influence on the total weight 
and, consequently, also no significant influence on the condition number. Actually, 
W* and W^ yielded precisely the same total weight for all functions that we have 
studied. Using that heuristics, the run time of Provan's algorithm improves to 

0(|E| + |V|log|V|), 

because its main part reduces to applying Dijkstra's shortest path algorithm just 
once. Fig. 2 compares W* and W^ for a few examples. 

3.3. Size of the Grid Domain. We restrict ourselves to graphs given by finite 
square grids of step size h, centered at z = — with all vertices and edges removed 
that do not belong to the domain D. Since Provan's algorithm just requires an 
embedded but not a planar graph, we may add the diagonals of the grid cells 
as further edges to the graph. These additional edges are advantageous, e.g., in 
all those cases which sensitively depend on properly matching the direction of 
steepest descent in the saddle points of d(z) (Bornemann 2011, §9). 

The side length / of the square domain (that is to be subdivided by the grid) has 
to be chosen large enough to contain a SEW that would approximate an optimal 
general integration contour. For example, if / is entire, we choose / > 2r*, where 
r* is the radius of the optimal circular contour given in Bornemann (2011). In 
other cases we employ a simple search for a suitable value of / by calculating W* 
for increasing values of I until d(W*) does not substantially decrease anymore. 
During this search each grid uses just a fixed number of vertices. 

3.4. Multilevel Refinement of the SEW. Choosing a proper value of h is not 
straightforward since we would like to balance a good approximation of a generally 
optimal integration contour with a reasonable amount of computing time. In 
principle, we construct a sequence of SEWs for smaller and smaller values of 
h until the weight of W* does not substantially decrease anymore. To avoid an 
unduly amount of computational work, we do not refine the grid everywhere 
but use an adaptive refinement by confining it to a tubular neighborhood of the 
currently given SEW W* (see Fig. 3): 

1: calculate W* within an initial grid; 

2: subdivide each rectangle adjacent to W* into 4 rectangles; 

3: remove all other rectangles; 

4: calculate W* in the newly created graph. 
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f(z) = e z , n = 300 



/(z) = e~ z ,n = 300 





/(z) = cos(z), n = 300 



/(z) = Ai(z) , n = 300 





/(z) = l/r(z) , n = 2006 



/(z) = (l-z) n/2 ,n = 10 





Figure 2. W* (red) vs. W^ (blue): the color coding shows the size 
of logrf(z); with red for large values and green for small values. The 
thin black lines are level curves of logrf(z); the smallest level shown is 
the threshold, below of which the edges of W* do not contribute to the 
first couple of significant digits of the total weight. The plots illustrate 
that W* and W^ differ just in a small region well below that threshold; 
consequently, both walks yield about the same condition number. 
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Figure 3. Multilevel refinement of W* (f(z) = 1/T(z), n = 2006) 



As long as the total weight of W* decreases substantially, steps 2 to 4 are repeated. 
It is even possible to optimize that process further by not subdividing rectangles 
that just contain vertices or edges of W* having weights below a certain threshold. 

3.5. Quadrature Rule for the Cauchy Integral. Finally, after calculation of the 
SEW r = W*, the Cauchy integral (1) has to be evaluated by some accurate 
numerical quadrature. We decompose T into maximally straight line segments, 
each of which can be a collection of many edges. On each of those line segments we 
employ Clenshaw-Curtis quadrature in Chebyshev-Lobatto points. Additionally 
we neglect segments with a weight smaller than 10 -24 times the maximum weight 
of an edge of T, since such segments will not contribute to the result within 
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Figure 4. Illustration of the spectral accuracy of piecewise Clenshaw- 
Curtis quadrature on SEW contours for a function with a branch-cut 
singularity. For larger n, we observe a significant improvement by adding 
diagonals to the grid. We get to machine precision for n = 10 and loose 
about two digits for n = 300. Note that for optimized circular contours 
the loss would have been about 6 digits for n = 10 and about 15 digits 
for n = 300 (cf. Bornemann 2011, Thm. 4.7). 

Table 1. Condition numbers for some f(z): r* are the optimal radii 
given in Bornemann (2011); W* was calculated in all cases (except the 
last one) on a 70 x 70-grid with I = 3r* (in the last case I was found 
by the method of §3.3). For 1/T(z), the specific orders of differentiation 
n = 2006 and n = 10935 are two of those very rare orders that give 
exceptionally large condition numbers (cf. Bornemann 2011, Table 5). 



/CO 



n K(W*,n) K(C r ^,n) 



cosz 
Ai(z) 



300 
300 
300 
300 



2.0 
1.1 
1.2 
1.6 



1.0 
1.0 
1.0 
1.2 



i/r(z) 


300 


2.2 


1.6 


i/r(z) 


2006 


7.8 • 10 4 


4.7 • 10 4 


i/r(z) 


10935 


1.6 • 10 5 


1.4 • 10 5 


(l-z) 11 / 2 


10 


in 


5.0 • 10 5 



machine precision. This way we not only get spectral accuracy but also, in many 
cases, less nodes than the vanilla version of trapezoidal sums on a circular contour 
would need: Fig. 4 shows an example with the order n = 300 of differentiation 
but accurate solutions using just about 200 nodes which is well below what the 
sampling condition would require for circular contours (Bornemann 2011, §2.1). 
Of course, trapezoidal sums would also benefit from some recursive device that 
helps to neglect those nodes which do not contribute to the numerical result. 



4- 



Numerical Results 



Table 1 displays condition numbers of SEWs W* on rectangular grids as com- 
pared to the optimal circles C r * for a couple of functions; Fig. 5 shows some of the 
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Figure 5. JA^ (admissible paths: blue = rectangular, magenta = rectan- 
gular with diagonals enabled) vs. C r ^ (cyan) for some examples of Table 1: 
the color coding shows the size of log^(z); with red for large values and 
green for small values. The thin black lines are level curves of logrf(z); 
the smallest level shown is the threshold, below of which the edges of 
JA^ do not contribute to the first significant digits of the total weight. 



corresponding contours. For entire / we observe that W*, like the optimal circle C r ^, 
automatically traverses the saddle points of d(z). It was shown in Bornemann 
(2011, Thm. 10.1) that, for such /, the major contribution of the condition number 
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comes from these saddle points and that circles are (asymptotically, as n — >> oo) 
paths of steepest decent. Since W* can cross a saddle point only in a horizonal, 
vertical, or (if enabled) diagonal direction, slightly larger condition numbers have 
to be expected. Nevertheless, the order of magnitude is precisely matched. This 
match holds in cases where circles give a best possible condition number of ap- 
proximately 1, as well as in cases with exceptionally large condition numbers, such 
as for f(z) = l/r(z) with the exceptional orders of differentiation n = 2006 and 
n = 10935 (cf. Bornemann 2011, §10.4). 

For some non-entire /, however, optimized circles can be far from optimal 
in general: Bornemann (2011, Thm. 4.7) shows that the optimized circle C r# for 
functions / from the Hardy space H 1 with boundary values in C k/0C yields a lower 
condition number bound of the form 

K(C u ,n) >cn k+a ; 

for instance, f(z) = (1 — z) n/2 gives K(C r ^,n) ~ 0.16059 • n 13/2 ; the principal 
branch of that / has a branch cut at (1, 00) and W* gives significantly better 
condition numbers than C r# by automatically following the cut. 

Conclusion. We observe that contours optimized in the class of rectangular paths 
are a flexible tool covering different classes of functions in a completely algorithmic 
fashion; no deep theory is needed to let the computation run (the theory is only 
required to explain large condition numbers if they cannot be avoided, such as for 
the entire function 1/T and certain orders n of differentiation). 
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